

<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
  "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">


<html xmlns="http://www.w3.org/1999/xhtml">
  <head>
    <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
    
    <title>scitools.convergencerate &mdash; SciTools 0.8.3 documentation</title>
    
    <link rel="stylesheet" href="../../_static/default.css" type="text/css" />
    <link rel="stylesheet" href="../../_static/pygments.css" type="text/css" />
    
    <script type="text/javascript">
      var DOCUMENTATION_OPTIONS = {
        URL_ROOT:    '../../',
        VERSION:     '0.8.3',
        COLLAPSE_INDEX: false,
        FILE_SUFFIX: '.html',
        HAS_SOURCE:  true
      };
    </script>
    <script type="text/javascript" src="../../_static/jquery.js"></script>
    <script type="text/javascript" src="../../_static/underscore.js"></script>
    <script type="text/javascript" src="../../_static/doctools.js"></script>
    <script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
    <link rel="top" title="SciTools 0.8.3 documentation" href="../../index.html" />
    <link rel="up" title="Module code" href="../index.html" /> 
  </head>
  <body>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="../../genindex.html" title="General Index"
             accesskey="I">index</a></li>
        <li class="right" >
          <a href="../../py-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li class="right" >
          <a href="../../np-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li><a href="../../index.html">SciTools 0.8.3 documentation</a> &raquo;</li>
          <li><a href="../index.html" accesskey="U">Module code</a> &raquo;</li> 
      </ul>
    </div>  

    <div class="document">
      <div class="documentwrapper">
        <div class="bodywrapper">
          <div class="body">
            
  <h1>Source code for scitools.convergencerate</h1><div class="highlight"><pre>
<span class="c">#!/usr/bin/env python</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">Module for estimating convergence rate of numerical algorithms,</span>
<span class="sd">based on data from a set of experiments.</span>

<span class="sd">Recommended usage:</span>
<span class="sd">Vary only discretization parameter (h in spatial problems, h/dt**q for</span>
<span class="sd">some q so that h/dt**q=const, in time-dependent problems with time step dt).</span>
<span class="sd">Use class OneDiscretizationPrm, or the function convergence_rate,</span>
<span class="sd">or the analyze_filedata convenience function. Start with reading</span>
<span class="sd">the convergence_rate function (too see an easily adapted example).</span>

<span class="sd">(There is support for fitting more general error models, like</span>
<span class="sd">C1*h**r1 + C2*h*dt**r2, with nonlinear least squares, but sound</span>
<span class="sd">fits are more challenging to obtain.)</span>
<span class="sd">&quot;&quot;&quot;</span>

<span class="kn">from</span> <span class="nn">Scientific.Functions.LeastSquares</span> <span class="kn">import</span> <span class="n">leastSquaresFit</span>
<span class="kn">from</span> <span class="nn">numpy</span> <span class="kn">import</span> <span class="n">zeros</span><span class="p">,</span> <span class="n">array</span><span class="p">,</span> <span class="n">asarray</span><span class="p">,</span> <span class="n">log10</span><span class="p">,</span> <span class="n">transpose</span><span class="p">,</span> <span class="n">linalg</span><span class="p">,</span> <span class="n">linspace</span>
<span class="kn">import</span> <span class="nn">sys</span>
<span class="kn">from</span> <span class="nn">scitools.std</span> <span class="kn">import</span> <span class="n">plot</span>

<span class="n">__all__</span> <span class="o">=</span> <span class="p">[</span><span class="s">&#39;ManyDiscretizationPrm&#39;</span><span class="p">,</span>
           <span class="s">&#39;OneDiscretizationPrm&#39;</span><span class="p">,</span>
           <span class="p">]</span>

<span class="n">log</span> <span class="o">=</span> <span class="n">log10</span>
<span class="n">inv_log</span> <span class="o">=</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="mi">10</span><span class="o">**</span><span class="n">x</span>

<span class="c"># The classes in this module have only static methods, i.e.,</span>
<span class="c"># classes are merely name spaces for related functions.</span>


<div class="viewcode-block" id="OneDiscretizationPrm"><a class="viewcode-back" href="../../convergencerate.html#scitools.convergencerate.OneDiscretizationPrm">[docs]</a><span class="k">class</span> <span class="nc">OneDiscretizationPrm</span><span class="p">(</span><span class="nb">object</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Tools for fitting an error model with only one discretization</span>
<span class="sd">    parameter: e = C*h^2.</span>
<span class="sd">    &quot;&quot;&quot;</span>

<div class="viewcode-block" id="OneDiscretizationPrm.error_model"><a class="viewcode-back" href="../../convergencerate.html#scitools.convergencerate.OneDiscretizationPrm.error_model">[docs]</a>    <span class="k">def</span> <span class="nf">error_model</span><span class="p">(</span><span class="n">p</span><span class="p">,</span> <span class="n">d</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;Return e = C*h**a, where p=[C, a] and h=d[0].&quot;&quot;&quot;</span>
        <span class="n">C</span><span class="p">,</span> <span class="n">a</span> <span class="o">=</span> <span class="n">p</span>
        <span class="n">h</span> <span class="o">=</span> <span class="n">d</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
        <span class="k">return</span> <span class="n">C</span><span class="o">*</span><span class="n">h</span><span class="o">**</span><span class="n">a</span>
</div>
<div class="viewcode-block" id="OneDiscretizationPrm.loglog_error_model"><a class="viewcode-back" href="../../convergencerate.html#scitools.convergencerate.OneDiscretizationPrm.loglog_error_model">[docs]</a>    <span class="k">def</span> <span class="nf">loglog_error_model</span><span class="p">(</span><span class="n">p</span><span class="p">,</span> <span class="n">d</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;As error_model if log-log data was used in the estimation.&quot;&quot;&quot;</span>
        <span class="n">C</span><span class="p">,</span> <span class="n">a</span> <span class="o">=</span> <span class="n">p</span>
        <span class="n">h</span> <span class="o">=</span> <span class="n">d</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
        <span class="k">return</span> <span class="n">log</span><span class="p">(</span><span class="n">C</span><span class="p">)</span> <span class="o">+</span> <span class="n">a</span><span class="o">*</span><span class="n">log</span><span class="p">(</span><span class="n">h</span><span class="p">)</span>

</div>
<div class="viewcode-block" id="OneDiscretizationPrm.linear_loglog_fit"><a class="viewcode-back" href="../../convergencerate.html#scitools.convergencerate.OneDiscretizationPrm.linear_loglog_fit">[docs]</a>    <span class="k">def</span> <span class="nf">linear_loglog_fit</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">e</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Linear least squares algorithm.</span>
<span class="sd">        Suitable for problems with only one distinct</span>
<span class="sd">        discretization parameter.</span>
<span class="sd">        *d* is the sequence of discretization parameter values, and</span>
<span class="sd">        *e* is the sequence of corresponding error values.</span>

<span class="sd">        The error model the data is supposed to fit reads</span>
<span class="sd">        ``log(e[i]) = log(C[i]) + a*log(d[i])``.</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="n">A</span> <span class="o">=</span> <span class="n">transpose</span><span class="p">(</span><span class="n">array</span><span class="p">([</span><span class="n">d</span><span class="p">,</span> <span class="n">zeros</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">d</span><span class="p">))</span><span class="o">+</span><span class="mi">1</span><span class="p">]))</span>
        <span class="n">sol</span> <span class="o">=</span> <span class="n">linalg</span><span class="o">.</span><span class="n">lstsq</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">e</span><span class="p">)</span>
        <span class="n">a</span><span class="p">,</span> <span class="n">logC</span> <span class="o">=</span> <span class="n">sol</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
        <span class="n">C</span> <span class="o">=</span> <span class="n">inv_log</span><span class="p">(</span><span class="n">logC</span><span class="p">)</span>
        <span class="k">return</span> <span class="n">a</span><span class="p">,</span> <span class="n">C</span>
</div>
<div class="viewcode-block" id="OneDiscretizationPrm.nonlinear_fit"><a class="viewcode-back" href="../../convergencerate.html#scitools.convergencerate.OneDiscretizationPrm.nonlinear_fit">[docs]</a>    <span class="k">def</span> <span class="nf">nonlinear_fit</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">e</span><span class="p">,</span> <span class="n">p0</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        ======== ===========================================================</span>
<span class="sd">        Name     Description</span>
<span class="sd">        ======== ===========================================================</span>
<span class="sd">        d        list of values of the (single) discretization</span>
<span class="sd">                 parameter in each experiment:</span>
<span class="sd">                 d[i] provides the values of the discretization,</span>
<span class="sd">                 parameter in experiement no. i.</span>
<span class="sd">        e        list of error values; e = (e_1, e_2, ...),</span>
<span class="sd">                 e[i] is the error associated with the parameters d[i]</span>

<span class="sd">        p0       starting values for the unknown parameters vector</span>
<span class="sd">        return   a, C; a is the exponent, C is the factor in front.</span>
<span class="sd">        ======== ===========================================================</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">d</span><span class="p">)</span> <span class="o">!=</span> <span class="nb">len</span><span class="p">(</span><span class="n">e</span><span class="p">):</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;d and e must have the same length&#39;</span><span class="p">)</span>
        <span class="k">if</span> <span class="ow">not</span> <span class="nb">isinstance</span><span class="p">(</span><span class="n">d</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="p">(</span><span class="nb">float</span><span class="p">,</span><span class="nb">int</span><span class="p">)):</span>
            <span class="k">raise</span> <span class="ne">TypeError</span><span class="p">(</span><span class="s">&#39;d must be an array of numbers, not </span><span class="si">%s</span><span class="s">&#39;</span> <span class="o">%</span> \
                            <span class="nb">str</span><span class="p">(</span><span class="nb">type</span><span class="p">(</span><span class="n">d</span><span class="p">[</span><span class="mi">0</span><span class="p">])))</span>
        <span class="c"># transform d and e to the data format required by</span>
        <span class="c"># the Scientific package:</span>
        <span class="n">data</span> <span class="o">=</span> <span class="p">[]</span>
        <span class="k">for</span> <span class="n">d_i</span><span class="p">,</span> <span class="n">e_i</span> <span class="ow">in</span> <span class="nb">zip</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">e</span><span class="p">):</span>
            <span class="n">data</span><span class="o">.</span><span class="n">append</span><span class="p">(((</span><span class="n">d_i</span><span class="p">,)</span> <span class="p">,</span> <span class="n">e_i</span><span class="p">))</span>  <span class="c"># recall (a,) conversion to tuple</span>
        <span class="n">sol</span> <span class="o">=</span> <span class="n">leastSquaresFit</span><span class="p">(</span><span class="n">OneDiscretizationPrm</span><span class="o">.</span><span class="n">error_model</span><span class="p">,</span> <span class="n">p0</span><span class="p">,</span> <span class="n">data</span><span class="p">)</span>
        <span class="n">C</span> <span class="o">=</span> <span class="n">sol</span><span class="p">[</span><span class="mi">0</span><span class="p">][</span><span class="mi">0</span><span class="p">]</span>
        <span class="n">a</span> <span class="o">=</span> <span class="n">sol</span><span class="p">[</span><span class="mi">0</span><span class="p">][</span><span class="mi">1</span><span class="p">]</span>
        <span class="k">return</span> <span class="n">a</span><span class="p">,</span> <span class="n">C</span>

</div>
<div class="viewcode-block" id="OneDiscretizationPrm.pairwise_rates"><a class="viewcode-back" href="../../convergencerate.html#scitools.convergencerate.OneDiscretizationPrm.pairwise_rates">[docs]</a>    <span class="k">def</span> <span class="nf">pairwise_rates</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">e</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Compare convergence rates, where each rate is based on</span>
<span class="sd">        a formula for two successive experiments.</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">d</span><span class="p">)</span> <span class="o">!=</span> <span class="nb">len</span><span class="p">(</span><span class="n">e</span><span class="p">):</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;d and e must have the same length&#39;</span><span class="p">)</span>

        <span class="n">rates</span> <span class="o">=</span> <span class="p">[]</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="nb">len</span><span class="p">(</span><span class="n">d</span><span class="p">)):</span>
            <span class="k">try</span><span class="p">:</span>
                <span class="n">rates</span><span class="o">.</span><span class="n">append</span><span class="p">(</span> <span class="n">log</span><span class="p">(</span><span class="n">e</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">/</span><span class="n">e</span><span class="p">[</span><span class="n">i</span><span class="p">])</span><span class="o">/</span><span class="n">log</span><span class="p">(</span><span class="nb">float</span><span class="p">(</span><span class="n">d</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">1</span><span class="p">])</span><span class="o">/</span><span class="n">d</span><span class="p">[</span><span class="n">i</span><span class="p">])</span> <span class="p">)</span>
            <span class="k">except</span> <span class="p">(</span><span class="ne">ZeroDivisionError</span><span class="p">,</span> <span class="ne">OverflowError</span><span class="p">):</span>
                <span class="n">rates</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="mi">0</span><span class="p">)</span>
        <span class="c"># estimate C from the last data point:</span>
        <span class="k">try</span><span class="p">:</span>
            <span class="n">C</span> <span class="o">=</span> <span class="n">e</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">/</span><span class="n">d</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span><span class="o">**</span><span class="n">rates</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
        <span class="k">except</span><span class="p">:</span>
            <span class="n">C</span> <span class="o">=</span> <span class="mi">0</span>
        <span class="k">return</span> <span class="n">rates</span><span class="p">,</span> <span class="n">C</span>

</div>
<div class="viewcode-block" id="OneDiscretizationPrm.analyze"><a class="viewcode-back" href="../../convergencerate.html#scitools.convergencerate.OneDiscretizationPrm.analyze">[docs]</a>    <span class="k">def</span> <span class="nf">analyze</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">e</span><span class="p">,</span> <span class="n">initial_guess</span><span class="o">=</span><span class="bp">None</span><span class="p">,</span>
                <span class="n">plot_title</span><span class="o">=</span><span class="s">&#39;&#39;</span><span class="p">,</span> <span class="n">filename</span><span class="o">=</span><span class="s">&#39;tmp.ps&#39;</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Run linear, nonlinear and successive rates models.</span>
<span class="sd">        d: list/array of discretization parameter.</span>
<span class="sd">        e: errors corresponding to d.</span>
<span class="sd">        Plot results for comparison of the three approaches.</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="c"># convert to NumPy arrays:</span>
        <span class="n">d</span> <span class="o">=</span> <span class="n">asarray</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="nb">float</span><span class="p">);</span>  <span class="n">e</span> <span class="o">=</span> <span class="n">asarray</span><span class="p">(</span><span class="n">e</span><span class="p">,</span> <span class="nb">float</span><span class="p">)</span>

        <span class="c"># linear least squares fit:</span>
        <span class="n">a1</span><span class="p">,</span> <span class="n">C1</span> <span class="o">=</span> <span class="n">OneDiscretizationPrm</span><span class="o">.</span><span class="n">linear_loglog_fit</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">d</span><span class="p">),</span> <span class="n">log</span><span class="p">(</span><span class="n">e</span><span class="p">))</span>
        <span class="k">print</span> <span class="s">&quot;linear LS fit: const=</span><span class="si">%e</span><span class="s"> rate=</span><span class="si">%.1f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">C1</span><span class="p">,</span> <span class="n">a1</span><span class="p">)</span>

        <span class="c"># nonlinear least squares fit (no log-log):</span>
        <span class="n">a2</span><span class="p">,</span> <span class="n">C2</span> <span class="o">=</span> <span class="n">OneDiscretizationPrm</span><span class="o">.</span><span class="n">nonlinear_fit</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">e</span><span class="p">,</span> <span class="n">initial_guess</span><span class="p">)</span>
        <span class="k">print</span> <span class="s">&quot;nonlinear LS fit: const=</span><span class="si">%e</span><span class="s"> rate=</span><span class="si">%.1f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">C2</span><span class="p">,</span> <span class="n">a2</span><span class="p">)</span>

        <span class="c"># pairwise estimate of the rate:</span>
        <span class="n">rates</span><span class="p">,</span> <span class="n">C3</span> <span class="o">=</span> <span class="n">OneDiscretizationPrm</span><span class="o">.</span><span class="n">pairwise_rates</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">e</span><span class="p">)</span>
        <span class="n">a3</span> <span class="o">=</span> <span class="n">rates</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>
        <span class="k">print</span> <span class="s">&quot;pairwise fit: const=</span><span class="si">%e</span><span class="s"> rate=</span><span class="si">%.1f</span><span class="s">&quot;</span> <span class="o">%</span> <span class="p">(</span><span class="n">C3</span><span class="p">,</span> <span class="n">a3</span><span class="p">)</span>
        <span class="k">print</span> <span class="s">&quot;all rates:&quot;</span><span class="p">,</span> <span class="n">rates</span>

        <span class="k">if</span> <span class="n">C1</span> <span class="o">&lt;</span> <span class="mi">0</span> <span class="ow">or</span> <span class="n">C2</span> <span class="o">&lt;</span> <span class="mi">0</span> <span class="ow">or</span> <span class="n">C3</span> <span class="o">&lt;</span> <span class="mi">0</span><span class="p">:</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;Some fits give negative const value! Cannot plot.&#39;</span><span class="p">)</span>
            <span class="k">return</span>

        <span class="c"># else log plot:</span>
        <span class="n">log_d1</span> <span class="o">=</span> <span class="n">linspace</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">d</span><span class="p">[</span><span class="mi">0</span><span class="p">]),</span> <span class="n">log</span><span class="p">(</span><span class="n">d</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]),</span> <span class="mi">2</span><span class="p">)</span>
        <span class="n">log_e1</span> <span class="o">=</span> <span class="n">log</span><span class="p">(</span><span class="n">C1</span><span class="p">)</span> <span class="o">+</span> <span class="n">a1</span><span class="o">*</span><span class="n">log_d1</span>
        <span class="n">log_e2</span> <span class="o">=</span> <span class="n">log</span><span class="p">(</span><span class="n">C2</span><span class="p">)</span> <span class="o">+</span> <span class="n">a2</span><span class="o">*</span><span class="n">log_d1</span>
        <span class="n">log_e3</span> <span class="o">=</span> <span class="n">log</span><span class="p">(</span><span class="n">C3</span><span class="p">)</span> <span class="o">+</span> <span class="n">a3</span><span class="o">*</span><span class="n">log_d1</span>
        <span class="n">plot</span><span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">d</span><span class="p">),</span> <span class="n">log</span><span class="p">(</span><span class="n">e</span><span class="p">),</span> <span class="s">&#39;yo&#39;</span><span class="p">,</span>
             <span class="n">log_d1</span><span class="p">,</span> <span class="n">log_e1</span><span class="p">,</span> <span class="s">&#39;r-&#39;</span><span class="p">,</span>
             <span class="n">log_d1</span><span class="p">,</span> <span class="n">log_e2</span><span class="p">,</span> <span class="s">&#39;b-&#39;</span><span class="p">,</span>
             <span class="n">log_d1</span><span class="p">,</span> <span class="n">log_e3</span><span class="p">,</span> <span class="s">&#39;g-&#39;</span><span class="p">,</span>
             <span class="n">legend</span><span class="o">=</span><span class="p">(</span><span class="s">&#39;data&#39;</span><span class="p">,</span>
                     <span class="s">&#39;linear LS log-log fit: </span><span class="si">%.1f</span><span class="s">*h^</span><span class="si">%.1f</span><span class="s">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">C1</span><span class="p">),</span> <span class="n">a1</span><span class="p">),</span>
                     <span class="s">&#39;nonlinear LS fit: </span><span class="si">%.1f</span><span class="s">*h^</span><span class="si">%.1f</span><span class="s">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">C2</span><span class="p">),</span> <span class="n">a2</span><span class="p">),</span>
                     <span class="s">&#39;successive rates, last two experiments: </span><span class="si">%.1f</span><span class="s">*h^</span><span class="si">%.1f</span><span class="s">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="n">log</span><span class="p">(</span><span class="n">C3</span><span class="p">),</span> <span class="n">a3</span><span class="p">)),</span>
             <span class="c">#axis=[log_d1[-1], log_d1[0], 1.3*log(e[-1]), 1.5*log(e[0])],</span>
             <span class="n">hardcopy</span><span class="o">=</span><span class="n">filename</span><span class="p">)</span>
</div>
    <span class="n">analyze</span> <span class="o">=</span> <span class="nb">staticmethod</span><span class="p">(</span><span class="n">analyze</span><span class="p">)</span>

    <span class="n">error_model</span> <span class="o">=</span> <span class="nb">staticmethod</span><span class="p">(</span><span class="n">error_model</span><span class="p">)</span>
    <span class="n">loglog_error_model</span> <span class="o">=</span> <span class="nb">staticmethod</span><span class="p">(</span><span class="n">loglog_error_model</span><span class="p">)</span>
    <span class="n">linear_loglog_fit</span> <span class="o">=</span> <span class="nb">staticmethod</span><span class="p">(</span><span class="n">linear_loglog_fit</span><span class="p">)</span>
    <span class="n">nonlinear_fit</span> <span class="o">=</span> <span class="nb">staticmethod</span><span class="p">(</span><span class="n">nonlinear_fit</span><span class="p">)</span>
    <span class="n">pairwise_rates</span> <span class="o">=</span> <span class="nb">staticmethod</span><span class="p">(</span><span class="n">pairwise_rates</span><span class="p">)</span>

<span class="c"># convenience function:</span></div>
<span class="k">def</span> <span class="nf">convergence_rate</span><span class="p">(</span><span class="n">discretization_prm</span><span class="p">,</span> <span class="n">error</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    Given two lists/arrays with discretization parameters and</span>
<span class="sd">    corresponding errors in a numerical method (element no. i</span>
<span class="sd">    in the two lists must correspond to each other), this</span>
<span class="sd">    function assumes an error formula of the form E=C*d^r,</span>
<span class="sd">    where E is the error and d is the discretization parameter.</span>
<span class="sd">    The function returns C and r.</span>

<span class="sd">    Method used: OneDiscretizationPrm.pairwise_rates is called</span>
<span class="sd">    and the final r value is used for return. A check that</span>
<span class="sd">    the rates are converging (the last three) is done.</span>
<span class="sd">    &quot;&quot;&quot;</span>
    <span class="n">rates</span><span class="p">,</span> <span class="n">C</span> <span class="o">=</span> <span class="n">OneDiscretizationPrm</span><span class="o">.</span><span class="n">pairwise_rates</span><span class="p">(</span><span class="n">discretization_prm</span><span class="p">,</span> <span class="n">error</span><span class="p">)</span>
    <span class="c"># check that there is no divergence at the end of</span>
    <span class="c"># the series of experiments</span>
    <span class="n">differences</span> <span class="o">=</span> <span class="p">[</span><span class="n">rates</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">-</span> <span class="n">rates</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">rates</span><span class="p">)</span><span class="o">-</span><span class="mi">1</span><span class="p">)]</span>
    <span class="c"># the differences between the rates should decrease, at least</span>
    <span class="c"># toward the end</span>
    <span class="k">try</span><span class="p">:</span>
        <span class="k">if</span> <span class="n">differences</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">&gt;</span> <span class="n">differences</span><span class="p">[</span><span class="o">-</span><span class="mi">2</span><span class="p">]:</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;The pairwise convergence rates do not &#39;</span>\
                             <span class="s">&#39;decrease toward the end:</span><span class="se">\n</span><span class="si">%s</span><span class="s">&#39;</span> <span class="o">%</span> \
                             <span class="nb">str</span><span class="p">(</span><span class="n">rates</span><span class="p">))</span>
    <span class="k">except</span> <span class="ne">IndexError</span><span class="p">:</span>
        <span class="k">pass</span>  <span class="c"># not enough data to check the differences list</span>
    <span class="k">return</span> <span class="n">C</span><span class="p">,</span> <span class="n">rates</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span>

<div class="viewcode-block" id="ManyDiscretizationPrm"><a class="viewcode-back" href="../../convergencerate.html#scitools.convergencerate.ManyDiscretizationPrm">[docs]</a><span class="k">class</span> <span class="nc">ManyDiscretizationPrm</span><span class="p">(</span><span class="nb">object</span><span class="p">):</span>
    <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">    General tool for fitting an error model containing an</span>
<span class="sd">    arbitrary number of discretization parameters.</span>
<span class="sd">    The error is a weighted sum of each discretization parameter</span>
<span class="sd">    raised to a real expoenent. The weights and exponents are</span>
<span class="sd">    the unknown parameters to be fitted by a nonlinear</span>
<span class="sd">    least squares procedure.</span>
<span class="sd">    &quot;&quot;&quot;</span>

<div class="viewcode-block" id="ManyDiscretizationPrm.error_model"><a class="viewcode-back" href="../../convergencerate.html#scitools.convergencerate.ManyDiscretizationPrm.error_model">[docs]</a>    <span class="k">def</span> <span class="nf">error_model</span><span class="p">(</span><span class="n">p</span><span class="p">,</span> <span class="n">d</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        Evaluate the theoretical error model (sum of C*h^r terms):</span>
<span class="sd">        sum_i p[2*i]*d[i]**p[2*i+1]</span>

<span class="sd">        ====== =======================================================</span>
<span class="sd">        Name   Description</span>
<span class="sd">        ====== =======================================================</span>
<span class="sd">        p      sequence ofvalues of  parameters (estimated)</span>
<span class="sd">        d      sequence of values of (known) discretization parameters</span>
<span class="sd">        return error evaluated</span>
<span class="sd">        ====== =======================================================</span>

<span class="sd">        Note that ``len(p)`` must be ``2*len(d)`` in this model since</span>
<span class="sd">        there are two parameters (constant and exponent) for each</span>
<span class="sd">        discretization parameter.</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">p</span><span class="p">)</span> <span class="o">!=</span> <span class="mi">2</span><span class="o">*</span><span class="nb">len</span><span class="p">(</span><span class="n">d</span><span class="p">):</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;len(p)=</span><span class="si">%d</span><span class="s"> != 2*len(d)=</span><span class="si">%d</span><span class="s">&#39;</span> <span class="o">%</span> <span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">p</span><span class="p">),</span><span class="mi">2</span><span class="o">*</span><span class="nb">len</span><span class="p">(</span><span class="n">d</span><span class="p">)))</span>
        <span class="nb">sum</span> <span class="o">=</span> <span class="mi">0</span>
        <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="nb">len</span><span class="p">(</span><span class="n">d</span><span class="p">)):</span>
            <span class="nb">sum</span> <span class="o">+=</span> <span class="n">p</span><span class="p">[</span><span class="mi">2</span><span class="o">*</span><span class="n">i</span><span class="p">]</span> <span class="o">*</span> <span class="n">d</span><span class="p">[</span><span class="n">i</span><span class="p">]</span><span class="o">**</span><span class="n">p</span><span class="p">[</span><span class="mi">2</span><span class="o">*</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span>
        <span class="k">return</span> <span class="nb">sum</span>
</div>
<div class="viewcode-block" id="ManyDiscretizationPrm.nonlinear_fit"><a class="viewcode-back" href="../../convergencerate.html#scitools.convergencerate.ManyDiscretizationPrm.nonlinear_fit">[docs]</a>    <span class="k">def</span> <span class="nf">nonlinear_fit</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">e</span><span class="p">,</span> <span class="n">initial_guess</span><span class="p">):</span>
        <span class="sd">&quot;&quot;&quot;</span>
<span class="sd">        @param d: list of values of the set of discretization</span>
<span class="sd">              parameters in each experiment:</span>
<span class="sd">              d = ((d_1,d_2,d_3),(d_1,d_2,d_3,),...);</span>
<span class="sd">              d[i] provides the values of the discretization</span>
<span class="sd">              parameters in experiement no. i.</span>
<span class="sd">        @param e: list of error values; e = (e_1, e_2, ...):</span>
<span class="sd">              e[i] is the error associated with the parameters</span>
<span class="sd">              d[i].</span>
<span class="sd">        @param initial_guess: the starting value for the unknown</span>
<span class="sd">        parameters vector.</span>
<span class="sd">        @return: list of fitted paramters.</span>
<span class="sd">        &quot;&quot;&quot;</span>
        <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">d</span><span class="p">)</span> <span class="o">!=</span> <span class="nb">len</span><span class="p">(</span><span class="n">e</span><span class="p">):</span>
            <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s">&#39;len(d) != len(e)&#39;</span><span class="p">)</span>
        <span class="c"># transform d and e to the data format required by</span>
        <span class="c"># the Scientific package:</span>
        <span class="n">data</span> <span class="o">=</span> <span class="p">[]</span>
        <span class="k">for</span> <span class="n">d_i</span><span class="p">,</span> <span class="n">e_i</span> <span class="ow">in</span> <span class="nb">zip</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">e</span><span class="p">):</span>
            <span class="k">if</span> <span class="nb">isinstance</span><span class="p">(</span><span class="n">d_i</span><span class="p">,</span> <span class="p">(</span><span class="nb">float</span><span class="p">,</span> <span class="nb">int</span><span class="p">)):</span>
                <span class="n">data</span><span class="o">.</span><span class="n">append</span><span class="p">(((</span><span class="n">d_i</span><span class="p">,),</span> <span class="n">e_i</span><span class="p">))</span>
            <span class="k">else</span><span class="p">:</span>  <span class="c"># d_i is tuple, list, array, NumArray, ...</span>
                <span class="n">data</span><span class="o">.</span><span class="n">append</span><span class="p">((</span><span class="n">d_i</span><span class="p">,</span> <span class="n">e_i</span><span class="p">))</span>
        <span class="n">sol</span> <span class="o">=</span> <span class="n">leastSquaresFit</span><span class="p">(</span><span class="n">ManyDiscretizationPrm</span><span class="o">.</span><span class="n">error_model</span><span class="p">,</span>
                              <span class="n">initial_guess</span><span class="p">,</span> <span class="n">data</span><span class="p">)</span>
        <span class="c"># return list of fitted parameters (p in error_model)</span>
        <span class="c"># (sol[1] is a measure of the quality of the fit)</span>
        <span class="k">return</span> <span class="n">sol</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
</div>
    <span class="n">error_model</span> <span class="o">=</span> <span class="nb">staticmethod</span><span class="p">(</span><span class="n">error_model</span><span class="p">)</span>
    <span class="n">nonlinear_fit</span> <span class="o">=</span> <span class="nb">staticmethod</span><span class="p">(</span><span class="n">nonlinear_fit</span><span class="p">)</span>

</div>
<span class="k">def</span> <span class="nf">_test1</span><span class="p">():</span>
    <span class="sd">&quot;&quot;&quot;Single discretization parameter test.&quot;&quot;&quot;</span>
    <span class="kn">import</span> <span class="nn">random</span>
    <span class="n">random</span><span class="o">.</span><span class="n">seed</span><span class="p">(</span><span class="mi">1234</span><span class="p">)</span>
    <span class="n">n</span> <span class="o">=</span> <span class="mi">7</span>
    <span class="n">h</span> <span class="o">=</span> <span class="mi">1</span>
    <span class="n">e</span> <span class="o">=</span> <span class="p">[];</span>  <span class="n">d</span> <span class="o">=</span> <span class="p">[]</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">7</span><span class="p">):</span>
        <span class="n">h</span> <span class="o">/=</span> <span class="mf">2.0</span>
        <span class="n">error</span> <span class="o">=</span> <span class="n">OneDiscretizationPrm</span><span class="o">.</span><span class="n">error_model</span><span class="p">((</span><span class="mi">4</span><span class="p">,</span><span class="mi">2</span><span class="p">),</span> <span class="p">(</span><span class="n">h</span><span class="p">,))</span>
        <span class="n">error</span> <span class="o">+=</span> <span class="n">random</span><span class="o">.</span><span class="n">gauss</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="mf">0.1</span><span class="o">*</span><span class="n">error</span><span class="p">)</span>  <span class="c"># perturb data</span>
        <span class="n">d</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">h</span><span class="p">)</span>
        <span class="n">e</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">error</span><span class="p">)</span>
    <span class="n">OneDiscretizationPrm</span><span class="o">.</span><span class="n">analyze</span><span class="p">(</span><span class="n">d</span><span class="p">,</span> <span class="n">e</span><span class="p">,</span> <span class="n">initial_guess</span><span class="o">=</span><span class="p">(</span><span class="mi">3</span><span class="p">,</span><span class="mi">2</span><span class="p">))</span>

<span class="k">def</span> <span class="nf">analyze_filedata</span><span class="p">():</span>
    <span class="c"># read filename and initial guess of C and r in error formula E=C*h^r</span>
    <span class="n">f</span> <span class="o">=</span> <span class="nb">open</span><span class="p">(</span><span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="s">&#39;r&#39;</span><span class="p">)</span>
    <span class="n">C</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">[</span><span class="mi">2</span><span class="p">])</span>
    <span class="n">r</span> <span class="o">=</span> <span class="nb">float</span><span class="p">(</span><span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">[</span><span class="mi">3</span><span class="p">])</span>
    <span class="k">try</span><span class="p">:</span>
        <span class="n">plot_title</span> <span class="o">=</span> <span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">[</span><span class="mi">4</span><span class="p">]</span>
    <span class="k">except</span><span class="p">:</span>
        <span class="n">plot_title</span> <span class="o">=</span> <span class="s">&#39;&#39;</span>
    <span class="kn">from</span> <span class="nn">scitools.filetable</span> <span class="kn">import</span> <span class="n">read</span>
    <span class="n">data</span> <span class="o">=</span> <span class="n">read</span><span class="p">(</span><span class="n">f</span><span class="p">)</span>
    <span class="k">print</span> <span class="n">data</span>
    <span class="n">OneDiscretizationPrm</span><span class="o">.</span><span class="n">analyze</span><span class="p">(</span><span class="n">data</span><span class="p">[:,</span><span class="mi">0</span><span class="p">],</span> <span class="n">data</span><span class="p">[:,</span><span class="mi">1</span><span class="p">],</span>
                                 <span class="n">initial_guess</span><span class="o">=</span><span class="p">(</span><span class="n">C</span><span class="p">,</span><span class="n">r</span><span class="p">),</span>
                                 <span class="n">plot_title</span><span class="o">=</span><span class="n">plot_title</span><span class="p">)</span>

<span class="c"># extensions:</span>
<span class="c"># example with dx, dy and dt</span>
<span class="c"># same example, but with factors to get a common rate</span>
<span class="c"># dx, dt tables and experiments with whole table, one</span>
<span class="c"># column and one row, and the diagonal</span>

<span class="k">if</span> <span class="n">__name__</span> <span class="o">==</span> <span class="s">&#39;__main__&#39;</span><span class="p">:</span>
    <span class="k">if</span> <span class="nb">len</span><span class="p">(</span><span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">)</span> <span class="o">==</span> <span class="mi">1</span><span class="p">:</span>
        <span class="k">print</span> <span class="s">&#39;Usage: </span><span class="si">%s</span><span class="s"> filename C-guess r-guess [plot-title]&#39;</span> <span class="o">%</span> <span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span>
    <span class="k">elif</span> <span class="n">sys</span><span class="o">.</span><span class="n">argv</span><span class="p">[</span><span class="mi">1</span><span class="p">]</span> <span class="o">==</span> <span class="s">&#39;example&#39;</span><span class="p">:</span>
        <span class="n">_test1</span><span class="p">()</span>
    <span class="k">else</span><span class="p">:</span>
        <span class="n">analyze_filedata</span><span class="p">()</span>
</pre></div>

          </div>
        </div>
      </div>
      <div class="sphinxsidebar">
        <div class="sphinxsidebarwrapper">
            <p class="logo"><a href="../../index.html">
              <img class="logo" src="../../_static/scitools_logo.jpg" alt="Logo"/>
            </a></p>
<div id="searchbox" style="display: none">
  <h3>Quick search</h3>
    <form class="search" action="../../search.html" method="get">
      <input type="text" name="q" />
      <input type="submit" value="Go" />
      <input type="hidden" name="check_keywords" value="yes" />
      <input type="hidden" name="area" value="default" />
    </form>
    <p class="searchtip" style="font-size: 90%">
    Enter search terms or a module, class or function name.
    </p>
</div>
<script type="text/javascript">$('#searchbox').show(0);</script>
        </div>
      </div>
      <div class="clearer"></div>
    </div>
    <div class="related">
      <h3>Navigation</h3>
      <ul>
        <li class="right" style="margin-right: 10px">
          <a href="../../genindex.html" title="General Index"
             >index</a></li>
        <li class="right" >
          <a href="../../py-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li class="right" >
          <a href="../../np-modindex.html" title="Python Module Index"
             >modules</a> |</li>
        <li><a href="../../index.html">SciTools 0.8.3 documentation</a> &raquo;</li>
          <li><a href="../index.html" >Module code</a> &raquo;</li> 
      </ul>
    </div>
    <div class="footer">
        &copy; Copyright 2012, H. P. Langtangen, J. Ring, ++.
      Created using <a href="http://sphinx.pocoo.org/">Sphinx</a> 1.1.2.
    </div>
  </body>
</html>